load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin

  plotname = "figure_NUMICE_acpd_v5_iop_v2" 
  plotform = "eps" 

  fna = "../diag_milled_acpd_v5.nc"
  fnb = "../diag_milled_acpd_v5_iop_v2.nc"
  fnc = "../sgp_20100401_shaocheng.nc"

  fla = addfile(fna,"r")
  flb = addfile(fnb,"r")

  xa = fla->NUMICE
  xb = flb->NUMICE
  xt = fla->T
  xp = fla->PS
  hyam = fla->hyam
  hybm = fla->hybm

  nt   = dimsizes(xa(:,0,0,0))
  nlev = dimsizes(xb(0,:,0,0))

  v1 = xa(:,:,0,0)
  v2 = xb(:,:,0,0)
  vt = xt(:,:,0,0)
  vpp = xp(:,0,0)

  rho = vt
  pre = vt

  do ik = 0,nlev-1
     pre(:,ik) = (/ doubletofloat(hyam(ik)) + doubletofloat(hybm(ik)) * vpp(:) /)
  end do

  r = 287.042 ;; J/kg/K.

  rho = pre / (r * vt)

  ;; #/kg -> #/L 

  v1 = v1 * rho / 1.e3 
  v2 = v2 * rho / 1.e3 

  vp = new((/2,nlev,nt/),"float")


 wks = gsn_open_wks(plotform,plotname)     ; open a ncgm file
 gsn_define_colormap(wks,"BlueYellowRed") ; choose colormap

; Individual plot resources
  
  res = True
  res@gsnDraw                = False
  res@gsnFrame               = False
  res@cnFillMode = "RasterFill"
 
  res@lbBoxMinorExtentF     = 0.15    
  res@lbOrientation         = "Vertical"
  res@lbLabelsOn            = True
  res@lbLabelFontHeightF    = 0.015 ;;2
  ;;res@lbLabelAutoStride     = True
  res@lbLabelStride  = 1
  res@lbTitleFontHeightF    = 0.015
  res@tiXAxisFontHeightF    = 0.02
  res@tiXAxisString         = "Days" 
  res@tiYAxisFontHeightF    = 0.02
  res@tiYAxisString         = "Pressure (mb)"
   
;  res@tmXBMinorOn           = False
  res@tmYLMinorOn           = False
  res@tmXTMinorOn           = False
  res@tmYRMinorOn           = False
  res@tmXTOn                = False
  
  ;;res@gsnRightString        = "" 
  ;;res@gsnLeftString         = "" 
  res@gsnStringFontHeightF  = 0.02
  res@gsnCenterStringFontHeightF  = 0.03
  
  res@cnFillOn             = True          ; turn on color
  res@gsnSpreadColors      = True     ; use full colormap
  res@cnLinesOn            = False    ; no contour lines
 
  res@gsnSpreadColors      = True          ; use full range of colormap
  
  ;;if (cmin.ne.cmax) then 
  ;;  res@cnLevelSelectionMode = "ManualLevels"
  ;;  res@cnLevelSpacingF      = cgap
  ;;  res@cnMinLevelValF       = cmin
  ;;  res@cnMaxLevelValF       = cmax
  ;;end if
  res@vpWidthF             = 0.8          ; change aspect ratio of plot
  res@vpHeightF            = 0.3
  ;;res@trYReverse           = True
  
  v2!0 = "time"  ;Reconstruct meta data as the time dimension has changed.
  v2!1 = "lev"
  
  v2&time = xa&time
  v2&lev  = xa&lev
  v2&lev@units = "hPa"
  v2&lev@long_name = "Pressure Levels"

  
  v1!0 = "time"
  v1!1 = "lev"
  
  v1&time = xa&time
  v1&lev  = xa&lev
  v1&lev@units = "hPa"
  v1&lev@long_name = "Pressure Levels"


  res@tmYRMode             = "Automatic"          ; turn off special labels on right axis

  plot = new(4,graphic)   


  res@cnLevelSelectionMode = "ExplicitLevels"    			; set explicit contour levels
  res@cnLevels    = (/ 1.e5,2.e5,5.e5,1.e6,2.e6,5.e6,1.e7,2.e7,5.e7,1.e8 /) / 1.e6 

  res@gsnLeftString = "A. IOP V1" 
  res@gsnRightString = "#/L" 

  plot(0) = gsn_csm_pres_hgt(wks,v1({lev|100:1000},time|:),res)

  res@cnLevels    = (/ 1.e5,2.e5,5.e5,1.e6,2.e6,5.e6,1.e7,2.e7,5.e7,1.e8 /) / 1.e6 
  res@gsnLeftString = "B. IOP V2" 
  res@gsnRightString = "#/L" 
  plot(1) = gsn_csm_pres_hgt(wks,v2({lev|100:1000},time|:),res)

  delete(res@cnLevels) 
  v3 = v1 
  v3 = v2 - v1 
  res@cnLevels    = (/ -5.e7,-2.e7,-1.e7,-5.e6,-2.e6,-1.e6,1.e6,2.e6,5.e6,1.e7,2.e7,5.e7 /) / 1.e6 

  res@gsnLeftString = "C. Abs. Diff" 
  res@gsnRightString = "#/L" 
  plot(2) = gsn_csm_pres_hgt(wks,v3({lev|100:1000},time|:),res)


  v4 = v1 
  v4 = (v2 - v1)/(v1+v2+1.e-20) * 100. 

  v4 = where((v1+v2).gt.1.,v4,-999) 
  v4@_FillValue = -999 

  delete(res@cnLevels) 

  res@cnLevels = (/ -200., -100.,-50.,-20.,-10.,-5., 0., 5.,10.,20.,50.,100.,200./) 

print("P3") 

  res@gsnLeftString = "Rel. Diff" 
  res@gsnRightString = "%" 
  plot(3) = gsn_csm_pres_hgt(wks,v4({lev|100:1000},time|:),res)


;; draw panel with white space added
 resP                 = True
 resP@gsnPanelYWhiteSpacePercent = 5
 resP@gsnPanelXWhiteSpacePercent = 5
 gsn_panel(wks,plot,(/3,1/),resP)


end







